Skyrmion Dynamics in Multiferroic Insulator 



(N 

o 

(N 

Cli 
CD 

c/3 



Ye-Hua Liu, 1 You-Quan Li, 1 and Jung Hoon Han 2 ' 3 '* 

1 Zhejiang Institute of Modern Physics and Department of Physics, 
Zhejiang University, Hangzhou 310027, People's Republic of China 
a Department of Physics and BK21 Physics Research Division, 
Sungkyunkwan University, Suwon 440-746, Korea 
3 Asia Pacific Center for Theoretical Physics, Pohang, Gyeongbuk 790-784, Korea 

(Dated: September 17, 2012) 

Recent discovery of Skyrmion crystal phase in insulating multiferroic compound Cu20Se03 calls 
for new ways and ideas to manipulate the Skyrmions in the absence of spin transfer torque from 
the conduction electrons. It is shown here that the position-dependent electric field, pointed along 
the direction of the average induced dipole moment of the Skyrmion, can induce the Hall motion 
of Skyrmion with its velocity orthogonal to the field gradient. Finite Gilbert damping produces 
longitudinal motion. We find a rich variety of resonance modes excited by a.c. electric field. 

PACS numbers: 75.85.+t, 75.70.Kw, 76.50.+g 
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Skyrmions are increasingly becoming commonplace 
sightings among spiral magnets including the metallic 
B20 compounds [1-5] and most recently, in a multiferroic 
insulator Cu20Se03[6]. Both species of compounds dis- 
play similar thickness-dependent phase diagrams[5, 6] de- 
spite their completely different electrical properties, high- 
lighting the generality of the Skyrmion phase in spiral 
magnets. Along with the ubiquity of Skyrmion matter 
comes the challenge of finding means to control and ma- 
nipulate them, in a device-oriented manner akin to efforts 
in spintronics community to control the domain wall and 
vortex motion by electrical current. Spin transfer torque 
(STT) is a powerful means to induce fast domain wall 
motion in metallic magnets [7, 8]. Indeed, current-driven 
Skyrmion rotation[9] and collective drift [10], originating 
from STT, have been demonstrated in the case of spiral 
magnets. Theory of current-induced Skyrmion dynam- 
ics has been worked out in Refs. [11, 12]. In insulating 
compounds such as Cu20Se03, however, the STT-driven 
mechanism does not work due to the lack of conduction 
electrons. 

As with other magnetically driven multiferroic 
compounds [13], spiral magnetic order in Cu20Se03 is 
accompanied by finite electric dipole moment. Recent 
work by Seki et al. [14] further confirmed the mecha- 
nism of electric dipole moment induction in Cu2 0Se03 
to be the so-called pd-hybridization[15-17]. In short, the 
pe?-hybridization mechanism claims the dipole moment 
Pjj for every oxygen-TM (transition metal) bond propor- 
tional to (Si ■ iij) 2 eij where i and j stand for TM and 
oxygen sites, respectively, and is the unit vector con- 
necting them. Carefully summing up the contributions 
of such terms over a unit cell consisting of many TM- 
O bonds, Seki et al. were able to deduce the dipole 
moment distribution associated with a given Skyrmionic 
spin configuration[14]. It is interesting to note that the 
numerical procedure performed by Seki et al. is pre- 
cisely the coarse-graining procedure which, in the text- 



book sense of statistical mechanics, is tantamount to the 
Ginzburg-Landau theory of order parameters. Indeed we 
can show that Seki et aVs result for the dipole moment 
distribution is faithfully reproduced by the assumption 
that the local dipole moment Pi is related to the local 
magnetization Si by 
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with some coupling A. A similar expression was pro- 
posed earlier in Refs. [18, 19] as the GL theory 
of Ba2CoGe2Oy[20], another known pd- hybridization- 
originated multiferroic material with cubic crystal struc- 
ture. Each site i corresponds to one cubic unit cell of 
Cu20Se03 with linear dimension a ~ 8.9A, and we have 
normalized Si to have unit magnitude. The dimension of 
the coupling constant is therefore [A] = C • m. 

Having obtained the proper coupling between dipole 
moment and the magnetizaiton vector in Cu20Se03 one 
can readily proceed to study the spin dynamics by solv- 
ing Landau-Lifshitz-Gilbert (LLG) equation. Very small 
values of Gilbert damping parameter are assumed in the 
simulation as we are dealing with an insulating magnet. 
A new, critical element in the simulation is the term aris- 
ing from the dipolar coupling 
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where we have used the magneto-electric coupling ex- 
pression in Eq. (1). In essence this is a field-dependent 
(voltage-dependent) magnetic anisotropy term. The to- 
tal Hamiltonian for spin is given by H = i?HDM + 
Hme, where i?HDM consists of the Heisenberg and the 
Dzyaloshinskii-Moriya (DM) exchange and a Zeeman 
field term. Earlier theoretical studies showed -Hhdm to 
stabilize the Skyrmion phase[l, 21-24]. 
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Two field orientations can be chosen independently in 
experiments performed on insulating magnets. First, the 
direction of magnetic field B determines the plane, or- 
thogonal to B, in which Skyrmions form. Second, the 
electric field E can be applied to couple to the induced 
dipole moment of the Skyrmion and used as a "knob" to 
move it around. Three field directions used in Rcf. [14] 
and the induced dipole moment in each case are classified 
as (I) B || [001], P = 0, (II) B || [110], P || [001], and 
(III) B || [111], P || [111]. One can rotate the spin axis 
appearing in Eq. (1) accordingly so that the ^-direction 
coincides with the magnetic field orientation in a given 
setup and the x-direction with the crystallographic [TlO]. 
In each of the cases listed above we obtain the magneto- 
electric coupling, after the rotation, 
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In cases (II) and (III) the E-field is chosen parallel to the 
induced dipole moment P, Ej = £jP, to maximize the 
effect of dipolar coupling. In case (I) where there is no 
net dipole moment for Skyrmions we chose E || [001] to 
arrive at a simple magneto-electric coupling form shown 
above. 

Suppose now that the E-field variation is sufficiently 
slow on the scale of the lattice constant a to allow the 
writing of the continuum energy, 
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It is assumed that all variables behave identically along 
the thickness direction, of length d. The "dipolar charge" 
density pnir) couples to the electric field E(r) in the same 
way as the conventional electric charge does to the poten- 
tial field in electromagnetism. The analogy is also useful 
in thinking about the Skyrmion dynamics under the spa- 
tially varying E-field as we will show. The continuum 
form of dipolar charge density in Eq. 4 is 
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Division by the unit cell area a 2 ensures that /9o(r) has 
the dimension of areal density. Values for the ferromag- 
netic case, Sf = 1, is subtracted in writing down the def- 
inition (5) in order to isolate the motion of the Skyrmion 
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FIG. 1: (color online) (a) Skyrmion configuration and (b)-(d) 
the corresponding distribution of dipolar charge density for 
three magnetic field orientations as in Ref. 14. (b) B || [001] 
(c) B || [110] (d) B || [111]. For each case, electric field is 
chosen as E || [001], E || [001] and E || [111], respectively. See 
text for the definition of dipolar charge density. As schemati- 
cally depicted in (a), the Skyrmion executes a Hall motion in 
response to electric field gradient. 



relative to the ferromagnetic background. Due to the 
subtraction, the dipolar charge is no longer equivalent to 
the dipole moment of the Skyrmion. The distribution of 
dipolar charge density for the Skyrmion spin configura- 
tion in the three cases are plotted in Fig. 1. In case (I) 
the total dipolar charge is zero. In cases (II) and (III) 
the net dipolar charges are both negative with the re- 
lation, Qd^/Qe) 1 ^ = n/3/2, where Qd, of order unity, 
is obtained by integrating pd(i") over the space of one 
Skyrmion and divide the result by the number of spins 
iVsk inside the Skyrmion. If the field variation is slow on 
the scale of the Skyrmion, then the point-particle limit is 
reached by writing pd(?) = QoNsk J2j $( r ~ r j) where Vj 
spans the Skyrmion positions, and identical charge Qu is 
assumed for all the Skyrmions. We arrive at the "poten- 
tial energy" of the collection of Skyrmion particles, 



H M v = -XQnN Sk -Y,E(r j ). 



(6) 



A force acting on the Skyrmion will be given as the gra- 
dient Fi = — Vii?ME- Inter-Skyrmion interaction is ig- 
nored. 

The response of Skyrmions to a given force, on the 
other hand, is that of an electric charge in strong 
magnetic field, embodied in the Berry phase action 
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(—2irSft,Qskd/a 3 )J2 : jf dt(rj x i\,) • z, where Qsk is the 
quantized Skyrmion charge[12, 25], and S is the size of 
spin. Equation of motion follows from the combination 
of the Berry phase action and Eq. (6), 



(7) 



where Vj is the j-th Skyrmion velocity. Typical Hall ve- 
locity can be estimated by replacing |V.E| with AE/lsu, 
where AE is the difference in the field strength between 
the left and the right edge of the Skyrmion and Isk is its 
diameter. Taking a 2 iVsk ~ 'gk we nnc ^ tne velocity 
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which gives the estimated drift velocity of 1 mm/s for the 
field strength difference of 10 3 V/m across the Skyrmion. 
Experimental input parameters of ?sk = 10~ 7 m, and 
A = 10 _32 C • m were taken from Rcf. [14] in arriving at 
the estimation, as well as the dipolar and the Skyrmion 
charges Qd ~ — 1 and Qsk = — 1. We may estimate the 
maximum allowed drift velocity by equating the dipolar 
energy difference XAE across the Skyrmion to the ex- 
change energy J, also corresponding to the formation en- 
ergy of one Skyrmion [24] . The maximum expected veloc- 
ity thus obtained is enormous, ~ 10 m/s for J ~ lmeV, 
implying that with the right engineering one can achieve 
rather high Hall velocity of the Skyrmion. In an encour- 
aging step forward, electric field control of the Skyrmion 
lattice orientation in the Cu2 0Se03 crystal was recently 
demonstrated [26] . 

Results of LLG simulation is discussed next. To start, 
a sinusoidal field configuration Ei — Eq sin (2%Xi/L x ) is 
imposed on a rectangular L x x L y simulation lattice with 
L x much larger than the Skyrmion size. In the absence 
of Gilbert damping, a single Skyrmion placed in such 
an environment moved along the "equi-potential line" in 
the y-direction as expected from the guiding-center dy- 
namics of Eq. (7). In cases (II) and (III) where the 
dipolar charges are nonzero the velocity of the Skyrmion 
drift is found to be proportional to their respective dipo- 
lar charges Qd as shown in Fig. 2. The drift velocity 
decreased continuously as we reduced the field gradient, 
obeying the relation (7) down to the zero velocity limit. 
The dipolar charge is zero in case (I), and indeed the 
Skyrmion remains stationary for sufficiently smooth E- 
field gradient. Even for this case, Skyrmions can move 
for field gradient modulations taking place on the length 
scale comparable to the Skyrmion radius, for the reason 
that the forces acting on the positive dipolar charge den- 
sity blobs (red in Fig. 1(a)) are not completely canceled 
by those on the negative dipolar charge density blobs 
(blue in Fig. 1(a)) for sufficiently rapid variations of the 



field strength gradient. A small but non-zero drift veloc- 
ity ensues, as shown in Fig. 2. Longitudinal motion along 
the field gradient begins to develop with finite Gilbert 
damping, driving the Skyrmion center to the position of 
lowest "potential energy" E(r). For the Skyrmion lat- 
tice, imposing a uniform field gradient across the whole 
lattice may be too demanding experimentally, unless the 
magnetic crystal is cut in the form of a narrow strip the 
width of which is comparable to a few Skyrmion radii. 
In this case we indeed observe the constant drift of the 
Skyrmions along the length of the strip in response to 
the field gradient across it. The drift speed is still pro- 
portional to the field gradient, but about an order of 
magnitude less than that of an isolated Skyrmion under 
the same field gradient. We observed the excitation of 
breathing modes of Skyrmions when subject to a field 
gradient, and speculate that such breathing mode may 
interfere with the drift motion as the Skyrmions become 
closed-packed. 
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FIG. 2: (color online) Skyrmion position versus time for cases 
(I) through (III) for sinusoidal electric field modulation (see 
text) with the Skyrmion center placed at the maximum field 
gradient position. The average Hall velocities (in arbitrary 
units) in cases (II) and (III) indicated in the figure are ap- 
proximately proportional to the respective dipolar charges, in 
agreement with Eq. (7). A small velocity remains in case 
(I) due to imperfect cancelation of forces across the dipolar 
charge profile. 

Several movie files are included in the supplementary 
information. Il.gif and III.gif give Skyrmion motion for 
Ei = Eq sin (2irxi/L x ) on L x x L y — 66 x 66 lattice for 
magneto-electric couplings (II) and (III) in Eq. (3). Ill- 
Gilbert. gif gives the same E-field as III.gif, with finite 
Gilbert damping a — 0.2. I. gif describes the case (I) 
where the average dipolar charge is zero, with a rapidly 
varying electric field Ei = Eoshi(2TTXi/X x ) and X x com- 
parable to the Skyrmion radius. The case of a narrow 
strip with the field gradient across is shown in strip.gif. 

Mochizuki's recent simulation [2 7] revealed that inter- 
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nal motion of Skyrmions can be excited with the uniform 
a.c. magnetic held. Some of his predictions were con- 
firmed by the recent microwave measurement [29]. Here 
we show that uniform a.c. electric field can also ex- 
cite several internal modes due to the magneto-electric 
coupling. Time-localized electric field pulse was applied 
in the LLG simulation and the temporal response %(£) 
was Fourier analyzed, where the response function x(t) 
refers to (1/2) Ei([5r(t)] 2 -[Sf(*)] 2 ), (1/2) Ei([5f(t)] 2 - 
[Sf(t)] 2 ), and (V3/2)Ei[<S?(t)] a for cases (I) through 
(III), respectively. (In Mochizuki's work, the response 
function was the component of total spin along the a.c. 
magnetic field direction.) 

In case (I), the uniform electric field perturbs the ini- 
tial cylindrical symmetry of Skyrmion spin profile so that 
J2i([Sf(t)] 2 — [Sf{t)} 2 ) becomes non-zero and the over- 
all shape becomes elliptical. The axes of the ellipse then 
rotates counter-clockwise about the Skyrmion center of 
mass as illustrated in supplementary figure, E-mode.gif. 
There are two additional modes of higher energies with 
broken cylindrical symmetry in case (I), labeled XI and 
X2 in Fig. 3 and included as Xl-mode.gif and X2- 
mode.gif in the supplementary. The rotational direction 
of the Xl-mode is the same as in E-mode, while it is the 
opposite for X2-mode. 

As in Ref. [27], we find sharply defined breathing 
modes in cases (II) and (III) at the appropriate resonance 
frequency oj, in fact the same frequency at which the a.c. 
magnetic field excites the breathing mode. The verti- 
cal dashed line in Fig. 3 indicates the common breath- 
ing mode frequency. Movie file Bl-mode.gif shows the 
breathing mode in case (III). Additional, higher energy 
B2-mode (B2-mode.gif) was found in cases (II) and (III), 
which is the radial mode with one node, whereas the Bl 
mode is nodeless. 

In addition to the two breathing modes, E-mode and 
the two X-modes are excited in case (II) as well due 
to the partly in-plane nature of the spin perturbation, 
-(\E{t)/2)J2 l ([SH t )? ~ [Sf(t)} 2 )- In contrast, case 
(III), where the perturbation -(y/3XE(t)/2) £ JS?(i)] 2 
is purely out-of-plane, one only finds the B-modes. As 
a result, case (I) and (III) have no common resonance 
modes or peaks, while case (II) has all the peaks (though 
XI and X2 peaks are small). Compared to the magnetic 
field-induced resonances, a richer variety of modes are 
excited by a.c. electric field. In particular, the E-mode 
has lower excitation energy than the B-mode and has a 
sharp resonance feature, which should make its detection 
a relatively straightforward task. Full analytic solution 
of the excited modes [28] will be given later. 

In summary, motivated by the recent discov- 
ery of magneto-electric material Cu20SeC>3 exhibiting 
Skyrmion lattice phase, we have outlined the theory of 
Skyrmion dynamics in such materials. Electric field gra- 
dient is identified as the source of Skyrmion Hall motion. 
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FIG. 3: (color online) (a) Absorption spectra for a.c. uniform 
magnetic field as in Mochizuki's work (reproduced here for 
comparison), (b) Absorption spectra for a.c. uniform electric 
field in cases (I) through (III). In case (I) where there is no 
net dipolar charge we find three low-energy modes E, XI, and 
X2. For case (III) where the dipolar charge is finite we find 
Bl and B2 radial modes excited. Case (II) exhibits all five 
modes. Detailed description of each mode is given in the text. 



Several resonant excitations by a.c. electric field are iden- 
tified. 
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